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RESUMEN 


El resumen sera traducido al espanol por los editores. We present our recent results from numerical 
simulations of a magnetized flow in the vicinity of a black hole in the context of the collapsar model for GRBs. 
The simulations show that after an initial transient, the flow settles into a complex convolution of several 
distinct, time-dependent flow components including an accretion torus, its corona and outflow, an inflow and 
an outflow in the polar funnel. We focus on studying the nature and connection between these components, in 
particular between the inflows and related outflows. We find that rotational and MHD effects launch, accelerate, 
and sustain the outflows. We also find that an outflow can be formed even when the collapsing envelope has 
initially a very weak magnetic field and a very small angular momentum. Our main conclusion is that even for 
a relatively weak initial magnetic field and a slow rotation, a gravitational collapse of a stellar envelope can 
lead to formation of a very strong and very fast jet. 

ABSTRACT 

We present our recent results from numerical simulations of a magnetized flow in the vicinity of a black hole 
in the context of the collapsar model for GRBs. The simulations show that after an initial transient, the flow 
settles into a complex convolution of several distinct, time-dependent flow components including an accretion 
torus, its corona and outflow, an inflow and an outflow in the polar funnel. We focus on studying the nature 
and connection between these components, in particular between the inflows and related outflows. We find 
that rotational and MHD effects launch, accelerate, and sustain the outflows. We also find that an outflow can 
be formed even when the collapsing envelope has initially a very weak magnetic field and a very small angular 
momentum. Our main conclusion is that even for a relatively weak initial magnetic field and a slow rotation, 
a gravitational collapse of a stellar envelope can lead to formation of a very strong and very fast jet. 
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1. INTRODUCTION 

Gamma-ray hurts (GRBs) are associated with 
the huge release of energy in a matter of seconds. 
The collapsar model is one of most promising sce¬ 
narios to explain these as well as other properties of 
GRBs (Woosley 1993; Paczyhski 1998; MacEadyen & 
Woosley 1999; Popham, Woosley & Fryer 1999;Mac- 
Fadyen, Woosley & Heger 2001; Proga et al. 2003). 
In this scenario, the collapsed iron core of a massive 
star accretes gas at a high rate (~ IMqS"^) pro¬ 
ducing a large neutrino flux, a powerful outflow, and 
a GRB. Many breakthroughs were made in study¬ 
ing GRBs. For example, the association of long du¬ 
ration GRBs with stellar collapse was firmly con¬ 
firmed (Hjorth et al. 2003, Stanek et al. 2003). 
Nevertheless, basic properties of the GRB central 
engine are quite uncertain. This is because the phys¬ 
ical conditions in the central engine are extreme and 
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complex (e.g., gravitational field is very strong, the 
temperature and the mass and energy densities are 
very high; large scale supersonic/relativistic flows 
and small scale turbulent flows are physically con¬ 
nected) . Additionally, magnetic fields are most likely 
very important in determining the properties of the 
central engine. 

Effects of magnetic fields in the context of 
the collapsar model have been studied by a few 
groups (e.g., Mizuno, Yamada, Koide, & Shibata 
2004; Mizuno, et al. 2004; De Villiers, Staff, & 
Guyed 2005). The main focus of these studies is 
on 2 and 3 dimensional general relativistic mag¬ 
netohydrodynamic (GR MHD) simulations of jets 
launched self-consistently from accretion disks orbit¬ 
ing Schwarzschild or Kerr black holes. These are 
very important studies as they capture a few of the 
key elements of the central engine and soon may in¬ 
clude more elements such as sophisticated equation 
of state and neutrino physics. 
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Here, we present and discuss results from 2.5- 
dimensional, magnetohydrodynamic (MHD) simula¬ 
tions of the collapsar model using pseudo-Newtonian 
potential. These simulations (see also Proga et al. 
2003) are an extension of the work of Proga & Begel- 
man (2003, hereafter PB03) who studied MHD ac¬ 
cretion flows onto a black hole (BH). In particular, 
the collapsar simulations include a realistic equation 
of state, photodisintegration of bound nuclei and 
cooling due to neutrino emission. The simulations 
presented here are also an extension of collapsar sim¬ 
ulations by MacFadyen & Woosley (1999), as they 
include very similar neutrino physics and initial con¬ 
ditions but are in the MHD instead hydrodynamical 
(HD) limit. 

2. MODELS 

Proga et al.’s simulation begins after the 1.7 M 0 
iron core of a 25 Mq presupernova star has collapsed 
and follows the ensuing accretion of the 7 Mq helium 
envelope onto the central black hole formed by the 
collapsed iron core. A spherically symmetric progen¬ 
itor model is assumed, but the symmetry is broken 
by the introduction of a small, latitude-dependent 
angular momentum and a weak split-monopole mag¬ 
netic field. For more details, see Proga et al. (2003). 
We refer to the model presented in Proga et al. as 
run A. We will also present results from model B 
which is a rerun of model A with the numerical res¬ 
olution increased by a factor of 2 in the latitudinal 
direction (see section 3.2). Run B has also a higher 
resolution in the radial direction compared to run A; 
the resolution increase is a function of radius and is 
of a factor of 19 at the smallest radii and of a fac¬ 
tor of 1.015 at the largest radii. Run B differs also 
from run A in the way we set up the initial condi¬ 
tions. However, despite these differences the gross 
properties of the solution during the later phase of 
the evolution are very similar for both runs. The 
details of run B are in Proga (2005; in preparation). 

3. RESULTS 
3.1. Fiducial run 

Eigure 1 shows time sequence of logarithmic den¬ 
sity (top) and toroidal magnetic field maps (bottom) 
overplotted with the direction of the poloidal velocity 
from run A during the early phase of the evolution. 
The figure illustrates how after a transient episode 
of infall, lasting about 0.13 s, the gas with I w 2Rsc 
starts to pile up outside the black hole and to form 
a thick torus bounded by a centrifugal barrier near 
the rotation axis. Soon after the torus forms (i.e., 
within a couple of orbital times at the inner edge). 


the magnetic field is amplified by shear and the mag- 
netorotational instability (MRI, e.g., Balbus & Hal- 
wey 1991). We note that the third possible mecha¬ 
nism to increase the magnetic field, compression is 
negligible here because the initial field as well as the 
infalling gas are radial (e.g., see the top right panel 
in Fig. 2). The fast growth of the magnetic fields 
is exemplihed by the growth of the toroidal field as 
shown in the bottom row of panels in Fig. 1. We 
have verified that most of the inner torus is unsta¬ 
ble to MRI, and that our simulations have enough 
resolution to resolve, albeit marginally, the fastest 
growing MRI mode (see also Section 3.2 for the pre¬ 
sentation of the results from run B). 

The magnetic field effects drive the time evolu¬ 
tion of the torus including the mass accretion onto 
a BH. Another important effect of magnetic fields is 
that the torus produces a magnetized corona and an 
outflow (e.g., the middle panels in Fig. 1). The pres¬ 
ence of the corona and outflow is essential to the evo¬ 
lution of the inner flow at all times and the entire flow 
close to the rotational axis during the latter phase of 
the evolution. The outflow very quickly becomes suf¬ 
ficiently strong to overcome supersonically infalling 
gas (the radial Mach number in the polar funnel near 
the inner radius is ~ 5) and makes its way outward, 
reaching the outer boundary at t = 0.25 s (see below 
and fig. 3). Due to limited computing time, run A 
was stopped at t = 0.28215 s, which corresponds 
to 6705 orbits of the flow near the inner boundary. 
We expect the accretion to continue much longer, 
roughly the collapse timescale of the Helium enve¬ 
lope (~ 10 s), as in MacFadyen & Woosley (1999) 

Figure 2 shows the time evolution of several gross 
properties run A: the mass accretion rate through 
the inner boundary (top left panel), total magnetic 
energy (top right panel), neutrino luminosity (bot¬ 
tom left panel) and radial Poynting and kinetic flux 
along the polar axis at r = 190 Rs (bottom right 
panel). Unless otherwise stated, all quantities in this 
paper are in cgs units. 

Initially, during a precollapse phase, Ma stays 
nearly constant at the level of 5x10^2 g s-h Dur¬ 
ing this phase the zero-1 gas inside the initial helium 
envelope is accreted. Around t = 0.13 s, Mq rises 
sharply as the gas from the initial helium envelope 
reaches the inner boundary. However, this gas has 
non-zero I and a rotational supported torus and its 
corona and outflow form (as illustrated in Fig. 1), 
causing a drop in Ma after it reaches a maximum of 
2 X 10^3 g s-3 at t = 0.145 s. 

The accretion rate reaches a minimum of 6 x 
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IqSi g g-1 at t « 0.182 s and then fluctuates with a 
clear long-term increase. This increase is caused by 
the contribution from gas with I < 2Rsc, which is 
directly accreted (without need to transport 1) from 
outside the main body of the torus (flow component 
C in Fig. 4 exemplifies a low I inflow outside the 
torus). The total mass and angular momentum ac¬ 
creted onto the BH during run A (0.3 s) are 0.1 M 0 
and 3 x 10^® g cm^ s“^, respectively. 

The top left panel of Fig. 2 shows the time evo¬ 
lution of the total magnetic energy (integrated over 
the entire computational domain). The late phase 
of the time evolution of the magnetic energy is char¬ 
acteristic of weakly magnetized rotating accretion 
flows. Most of the magnetic energy is due to the 
toroidal component of the field. We note a huge in¬ 
crease of the toroidal magnetic field coinciding with 
the formation and development of the torus. Both 

and Bg are practically zero while B^ is constant 
during the precollapse phase of the evolution. But 
at t = 0.14 s the total energy in B^jy equals that in 
Br and just 0.025 s later the energy is higher 
than the Br energy by a factor of 50. At the end 
of simulations the total kinetic energy from the ra¬ 
dial, latitudinal and rotational motion are 4 x 10®®, 
6.5 X 10^®, and 2.3 x 10®^ erg, respectively. These 
gross properties indicate that the magnetic energy is 
large enough to play an important role in the flow 
dynamics. 

The bottom left panel of Fig. 2 shows the time 
evolution of the neutrino luminosity, Ly,. The neu¬ 
trino luminosity was computed under the assump¬ 
tion that all the gas in the model is optically thin 
to neutrinos. The neutrino emission stays at a rel¬ 
atively constant level of 3 x 10®® erg s“^ after the 
torus forms. This indicates that the gross properties 

- such as the temperature, density, and total mass 

- of the densest and hottest parts of the flow (the 
torus) also stay constant during the late phase of 
the evolution. 

The bottom right panel in Fig. 2 shows the area- 
integrated radial fluxes of magnetic and kinetic en¬ 
ergy at r = 190 Rs inside the polar outflow. For¬ 
mally, the polar outflow is defined as the region 
where > 0 and (3 < 1. The Poynting fluxes, stays 
nearly constant at the level of ~ 1 x 10®^ egr s“^ 
between t = 0.185 s (the time when the jet reached 
the radius of 190 Rs) and t = 0.22 s then fluctuates 
with a clear long-term decrease. Comparing the two 
fluxes, we conclude that the outflow is Poynting flux- 
dominated, with the Poynting flux exceeding the ki¬ 
netic energy flux by up to an order of magnitude. We 
note that the radial fluxes continue to vary with time 


by a factor up to 10 even during the late phase of 
the evolution and they are anti-correlated with the 
mass accretion. (The radial fluxes show much less 
detail compared to other properties shown in Fig. 2 
because there were computed less frequently during 
the course of the simulation.) 

The decrease of the radial fluxes and their anti¬ 
correlation with Ma is quite a surprising result be¬ 
cause the jet is supposed to be powered by accre¬ 
tion! However, the situation in a collapsing star is 
complex because Ma has two sources: (i) an accre¬ 
tion disk which does power a strong jet and (ii) the 
low I gas that accretes directly onto a BH and can 
prevent development of the torus jet. This infalling 
gas can produce an outflow/jet of its own but with 
properties different from those of the torus jet (see 
below). 

Figure 3 shows maps of the flow pattern and 
maps of Byf, on four different length scales at < = 
0.285 s (the left panels show the inner most part of 
the flow whereas the right panels show the flow on 
the largest scale). The main outflow is magnetically 
driven from the torus. Soon after the torus forms, 
the magnetic field very quickly deviates from its ini¬ 
tial radial configuration due to MRI and shear (see 
Fig. 2). This leads to fast growth of the toroidal 
magnetic field as field lines wind up due to the dif¬ 
ferential rotation. As a result the toroidal held dom¬ 
inates over the poloidal held and the gradient of the 
former drives an outflow. 

A comparison between the density and B^ maps 
shows that the polar regions of low density and high 
B^ coincide with the region of an outflow. Proga 
et al. (2003) noted that during the latter phase of 
the evolution not all of the material in the outflow 
originated in the innermost part of the torus - a 
significant part of the outflow is “peeled off” the in¬ 
falling gas at large radii by the magnetic pressure. 
However, as our later analysis and simulations show 
the outflow from the infalling gas is caused not only 
by the peeling off effect but it is also produced by the 
magnetic and thermal effects operating inside the in¬ 
falling gas itself. Thus there are not two but three 
types of outflows from a collapsing star (see Fig. 4). 

The first two panels of the left hand side of Fig. 3 
illustrate that the inner torus and its corona and out¬ 
flow cannot always prevent the low-1 gas from reach¬ 
ing the BH. Even the magnetic field cannot do it if 
the density of the incoming gas is too high (compare 
the upper and lower halfs of the panels in Figs 3 and 
5). 

We finish the presentation of the results from 
run A with an analysis of the properties of the ac- 
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cretion torus, i.e., the main engine in the model. 
Figure 4 shows the radial profiles of several quan¬ 
tities in run A, angle-averaged over a small wedge 
near the equator (between 6 = 86° and 94°), and 
time-averaged over 50 data files covering a period at 
the end of the simulations (from 0.2629 s through 
0.2818 s). We indicate the location of the last sta¬ 
ble circular orbit by the vertical dotted line in each 
panel. 

We measure the Reynolds stress, Ugas =< 
pVrSv^ > /P, and the Maxwell stress normalized 
to the magnetic pressure, amag =< 2BrB^/>. 
Note that Fig. 4 shows only the magnitude, not the 
sign, of the normalized stresses. We find that except 
for r ^ 2.5Rs and 10i?s ^ r ^ 12i?s, the Maxwell 
stress dominates over the Reynolds stress in the in¬ 
ner flow. The last panel in Fig. 4 shows that the 
toroidal component of the magnetic field is dominant 
for r < 50 i?s. 

Proga et al.’s simulation explored a relatively 
conservative case as they allowed for neutrino emis¬ 
sion but did not allow for the emitted neutrinos to 
interact with the gas or annihilate. The only sources 
of nonadiabatic heating in their simulation are the 
artificial viscosity and resistivity. 

3.2. High resolution run 

To explore the nature of the multicomponent 
flow resulting from a collapse of a rotating star we 
rerun model A using higher numerical resolution 
[see Section 3.1 and Proga (2005 in preparation)]. 
Higher resolution simulations are especially required 
to study the evolution and effects of MRI. Run B 
is qualitatively consistent with the results from the 
fiducial model. In particular, all the flow compo¬ 
nents found in run A can be also found in run B. 
Additionally, the properties of the flow components 
found in run A are consistent with those in run A. 
For example, MRI drives the evolution of the equa¬ 
torial torus and the simulations resolve the fastest 
growing MRI mode. 

Figure 5 shows the inner most part of the flow in 
run B at t=0.582 s. The seven main flow components 
are marked with upper case letters (see the figure 
caption). 

In the context of GRBs, the most important com¬ 
ponent of the flow is an outflow. The biggest differ¬ 
ence between the MHD collapsar model and the HD 
collapsar model considered by MacFadyen Woosley 
(1999) is that in the MHD limit, the outflow that 
develops soon after a torus forms is so strong that it 
breaks through a star even when very little of a low 
I gas is accreted whereas in the HD limit, a torus 


wind is relatively weak and the polar funnel must be 
evacuated if the wind were to leave the star. 

3.3. An outflow from a very I infall 

Generally, runs A and B show that large-scale 
magnetic fields can produce two types of outflows: 
(1) a jet from a rotationally supported accretion disk 
or torus and (2) an outflow from extremely low an¬ 
gular momentum gas that almost radially accretes 
onto a BH (e.g., PB03, Proga et al. 2003). Both an¬ 
alytic and numerical studies support the view that 
the former is robust, (e.g., Blandford & Payne 1982; 
Blandford 1990; De Villier et al. 2003). In fact, most 
MHD simulations of an accretion disk or torus show 
outflows (e.g., Uchida & Shibata 1985; Stone & Nor¬ 
man 1994; Hawley & Balbus 2002; De Villier et al. 
2004; Mizuno et al. 2004; Kato et al. 2004; McK¬ 
inney & Gammie 2004). However, the latter which 
is a simple, self-consistent solution for the MHD jet 
problem, has not been studied much. To articulate 
the basic physics that occurs in jet production, Proga 
(2005) performed simulations of a flow with angular 
momentum so low that, if not for the effects of MHD, 
the flow would accrete directly onto a BH without 
forming a disk. These simulations used simplified 
physics (i.e., no neutrino cooling and an adiabatic 
equation of state) similar to that explored by PB03. 

In Proga (2005), we found that even with a very 
weak initial magnetic field, the flow settles into a 
configuration with four components: (i) an equato¬ 
rial inflow, (ii) a bipolar outflow, (iii) a polar funnel 
outflow, and (iv) a polar funnel inflow. The sec¬ 
ond flow component of the MHD flow represents a 
simple yet robust example of a well-organized in¬ 
flow/outflow solution to the problem of MHD jet 
formation (see Fig. 6) and is the same in nature as 
components G and B in run B as well as in run A 
(e.g.. Fig. 5). The outflow from the low I infall is 
heavy, highly magnetized, and driven by magnetic 
and centrifugal forces. A significant fraction of the 
total energy in the jet is carried out by a large scale 
magnetic field. The properties of this outflow help 
to understand that the time evolution of the out¬ 
flow in run A. As we mentioned in section 3.1, the 
low I inflow can block or prevent development of the 
torus jet, the latter being light. The low I inflow 
produces the outflow that is heavy and although a 
significant fraction of its total energy is carried out 
by a large scale magnetic field it is not Poynting flux- 
dominated, contrary to the torus jet. 

A comparison between various simulations in¬ 
cluding those presented in section 3.1 and 3.2, where 
specific angular momentum was higher than that as- 
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sumed in Proga (2005), indicates that the flow com¬ 
ponents B and C develop for a wide range of the 
properties of the flow near the equator and near the 
poles. 

4. CONCLUSIONS 

Fully 3 dimensional GR MHD simulations are re¬ 
quired to capture many of the effects and instabili¬ 
ties of a magnetized fluid in a rotating star collapsing 
onto a BH . Neutrino physics, a sophisticated equa¬ 
tion of state and self-gravity should also be included. 
However, such full treatment of the MHD collapsar 
model is beyond reach of the current numerical codes 
at least for now. 

Here we present results from time-dependent 
two-dimensional MHD simulations of the collapsar 
model using pseudo-Newtonian potential. The simu¬ 
lations show that: 1) soon after the rotationally sup¬ 
ported torus forms, the magnetic field very quickly 
starts deviating from purely radial due to MRI and 
shear. This leads to fast growth of the toroidal 
magnetic field as field lines wind up due to the 
torus rotation; 2) The toroidal held dominates over 
the poloidal field and the gradient of the former 
drives a torus outflow against supersonically accret¬ 
ing gas through the polar funnel; 3) The torus out¬ 
flow is Poynting flux-dominated; 4) The torus out¬ 
flow reaches the outer boundary of the computa¬ 
tional domain (5 x 10® cm) with an expansion ve¬ 
locity of 0.2 c; 5) The torus outflow is in a form of 
a relatively narrow jet (when the jet breaks through 
the outer boundary its half opening angle is 5°); 6) 
Most of the energy released during the accretion is 
in neutrinos, = 2 x 10®^ erg s“^. Neutrino driv¬ 
ing will increase the outflow energy (e.g.. Fryer & 
Meszaros 2003 and references therein), but could also 
increase the mass loading of the outflow if the energy 
is deposited in the torus. A comparison of the MHD 
simulations with their HD counterparts show that 
a strong outflow breaks through a magnetized star 
much sooner than through a non-magnetized star. 

The above conclusions were reached by Proga et 
al. (2003) and here we confirmed them using a higher 
resolution simulation. However, we emphasize the 
fact that the flow settles into a complex convolution 
of several distinct, time-dependent flow components 
and the above mentioned torus and its outflow are 
just two of them. Other flow components include 
a torus corona and low I flows. A rotationally sup¬ 
ported torus and its corona and outflows were exten¬ 
sively studied in the past and remain a focus of many 
studies. We stress that in the context of the collap¬ 
sar model, where very low I and high I fluids are 


present, the situation is more complex. Therefore, 
future work, even without neutrino physics or effects 
of general relativity, is important to explore connec¬ 
tion and interaction of all the flow components, and 
their observational implications. 
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Fig. 1. Sequence of logarithmic density (top) and toroidal magnetic field maps (bottom) overplotted with the direction 
of the poloidal velocity from run A at times 0.153, 0.161, 0.165, and 0.173 s. The sequence illustrates the early phase 
of the formation of a rotational supported accretion torus and and of a magnetically driven outflow. The length scale 
is in units of the BH radius (i.e., r' = r/Rs and z' = z/Rs)- 
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Fig. 2. The time evolution of the mass accretion rate (top left panel), total magnetic energy due to each of the three 
field components (top right panel), neutrino luminosity (bottom left panel) and area-integrated radial Poynting and 
kinetic flux in the polar outflow at r = 190 Rs (bottom right panel) for run A. Formally, we define the polar outflow as 
the region where Ur > 0 and /3 < 1. Note the difference in the time range in the panel with the radial fluxes. 
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Fig. 3. Maps of logarithmic density (top) and toroidal magnetic field maps (bottom) overplotted with the direction of 
the poloidal velocity from rnn A at the end of simulations, i.e. time 0.2815 s. The length range increases from to the 
right. Note that the accretion torus is relatively small (it spans from 1 to about 20 Rs). Nevertheless, this tiny torus 
generates an outflow and mass and energy that can the dynamics and structure of the collapsing star over a large range 
of radii along the rotational axis. 
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Fig. 4. Radial profiles of various quantities from our run, time-averaged from 0.2629 through 0.2818 s. To construct 
each plot, we averaged the profiles over angle between Q = 86° and 94°. The top left panel plots the density (solid 
line) and temperature (dashed line). The top middle panel plots the gas pressure (solid line) and magnetic pressure. 
The top right panel plots the rotational, radial, Keplerian, and Alfven velocities (solid, dashed, dot-dashed, and dotted 
line, respectively), as well as the sound speed (triple-dot dashed line). The bottom left panel plots the angular velocity 
in units of 2c/i?s. The bottom middle panel plots the Maxwell stress, OLmag, and the Reynolds stress, a.gaa (solid and 
dashed line, respectively). We calculate the Reynolds stress using eq. (15) in PB03 and show only its amplitude. The 
bottom right panel plots the radial, latitudinal and toroidal components of the magnetic field (dot-dashed, dashed, and 
solid line, respectively). The length scale is in units of the BH radius (i.e., r' = rjRs). 
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Fig. 5. A map of logarithmic density overplotted with the direction of the poloidal velocity (left panel) and a contour 
map of specific angular momentum, I (right panel) from high resolution run B at t = 0.0582 s. The specific angular 
momentum is in units of 2Rsc. The minimum of I (contour closest to the rotational z-axis) is 0.2, and the contour 
levels are equally spaced at intervals of Z = 0.2. The maximum of I is 1.0 and its contour is plotted using dotted curves, 
whereas all the other contours are plotted using solid curves. This figure shows an inner most part of the flow when a 
torus just formed and started to develop an outflow. Note that the latter pushed aside the polar funnel accretion flow 
only below the equator. This figure was chosen to illustrate the complexity of the inner most part of the flow inside a 
collapsing magnetized star. The upper case letters, in the right panel, mark seven major components of this complex 
flow: A - a very low I polar funnel accretion flow; B - a highly magnetized outflow generated from a low I inflow C; C - 
an inflow with angular momentum so low that the flow would accrete directly onto a BH; D - a rotationally supported, 
MHD, turbulent accretion torus; E - a magnetized torus corona; F - a highly magnetized outflow generated from the 
torus; and G - a polar funnel outflow driven thermally (compare with A). See Proga (2005; in preparation) for more 
details. 
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Fig. 6. Maps of logarithmic density (left panel) and toroidal magnetic field (right panel) overplotted with an example 
of a streamline corresponding to an inflow/outflow for model C. The maps are also overplotted with the direction of the 
poloidal velocity and the direction of the poloidal field (the left and right panels, respectively). See Proga (2005) for 
more details. 
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